Quantum initial value representations using approximate Bohmian trajectories 
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Quantum trajectories, originating from the de Broglie-Bohm (dBB) hydrodynamic description 
of quantum mechanics, are used to construct time-correlation functions in an initial value rep- 
resentation (IVR). The formulation is fully quantum mechanical and the resulting equations for 
the correlation functions are similar in form to their semi-classical analogs but do not require the 
computation of the stability or monodromy matrix or conjugate points. We then move to a local 
trajectory description by evolving the cumulants of the wave function along each individual path. 
The resulting equations of motion are an infinite hierarchy, which we truncate at a given order. We 
show that time-correlation functions computed using these approximate quantum trajectories can 
be used to accurately compute the eigenvalue spectrum for various potential systems. 



I. INTRODUCTION 

Over the past few years there has been significant ef- 
fort in the development of the so-called hydrodynamic or 
quantum traiectory based approach to quantum mechan- 

ics. laaaiiiiiiaiEaiiiiiiEiiiiiiiiiiiia 

This formulation is based upon Bohm's "hidden variable" 
representation 0, 0, |23| is initiated by substituting the 
amplitude-phase decomposition of the time-dependent 
quantum wavefunction, ipiy^t) — R{y,t) exp{iS{y,t)/h), 
into the time-dependent Schrodinger equation. 
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Substituting Eq. 1 into Eq. 2 and separating the real and 
imaginary components yields equations of motion for p 
and S, 
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where p — is the probability density and S is the 
phase. These last two equations can be easily identi- 
fied as the continuity equation (Eq. ^ and the quan- 
tum Hamilton- Jacobi equation (Eq. In Eq. |31we can 
identify three contributions to the action. The first two 
are the kinetic and potential energies and the third is 
the quantum potential Q 0, . This term is best de- 
scribed as a shape kinetic energy since it is determined 
by the local curvature of the quantum density. In nu- 
merical applications, computation of the quantum po- 
tential is frequently rendered more accurate if deriva- 
tives are evaluated using the log-amplitude C = \og{R) 
[C is also referred to as the C-amplitude). In terms of 
derivatives of this amplitude, the quantum potential is 
Q = -?i2(v2c+ (vC)2)/2m. 

Eq. 13 and 121 are written in terms of partial derivatives 
in time. In order to compute the time-evolution of p (or 
equivalently C) and S along a give path, x{t) we need to 



move to a moving Lagrangian frame by transforming dt — 
dt + x(t)dx- The Lagrangian form of the hydrodynamic 
equations of motion resulting from this analysis are given 
by: 
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in which the derivative on the left side is appropriate 
for calculating the rate of change in a function along a 
fluid trajectory. Eq. |H| is a Newtonian-type equation in 
which the flow acceleration is produced by the sum of 
the classical force, fc — ~WV, and the quantum force is 
/, = -VQ 

The quantum trajectories obey two important non- 
crossing rules: (1) they cannot cross nodal surfaces (along 
which ■0 = 0); (2) they cannot cross each other. (In 
practice, because of numerical inaccuracies, these condi- 
tions may be violated.) Because of these two conditions, 
the quantum trajectories are very different from classical 
paths and represent a geometric optical rendering of the 
evolution of the quantum wave function. 

The primary difflculty in implementing a quantum tra- 
jectory based approach is in constructing the quantum 
potential. In spite of 3-4 years of intense effort by a 
number of groups 0, S S II II II 0, 111 111 111 El in- 
cluding our own, there has yet to emerge a satisfactory 
and robust way to compute Q and the quantum force 
— VQ directly from the quantum trajectories in multidi- 
mensional systems. This severely limits the applicability 
of a quantum trajectory based approach for computing 
exact quantum dynamics for all but the most trivial sys- 
tems. 

In this paper, we show that quantum trajectories can 
be used to construct exact quantum initial value rep- 
resentations (IVR) of correlation functions. The basic 
equations are extremely straightforward to derive and fol- 
low directly from a Bohm trajectory based description. 
The heart and sole of our approach, however, is the use 
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of an cummulant /derivative propagation scheme '2l'| that 
generates local equations of motion for a given quantum 
trajectory and all the derivatives of C and S along the 
path, allowing us to compute the correlation function on 
a trajectory by trajectory basis. In our numerical imple- 
mentation of this approach, we examine a highly anhar- 
monic one dimensional system (an inverted gaussian well) 
and find that the Bohni-IVR approach is surprisingly ac- 
curate in computing correlation functions even though 
at long-times the quantum trajectories themselves show 
crossings. 

II. THE BOHM INITIAL VALUE 
REPRESENTATION 



is the time-correlation function for some quantum oper- 
ator A following initial state tpo evolving under Hamilto- 
nian, H. 

Within the Bohmian framework, the synthesis of a 
quantum wavefunction along a quantum trajectory, x{t), 
is given by 

ijix,t) = J dx{t)6{x ~ x{t))J-^/^{x{t),x{0)) 

X e*^(^(°)'^(*»/''Vo(a;(0)). (8) 

Here, 



While in principle, the quantum wavefunction or the 
density matrix tells us everything we can know about a 
given physical process, often what we are interested in 
is a correlation function or an expectation value. This 
can be equivalently cast as either a time-correlation or a 
spectral function via Fourier transform: 

(Aiu;))=:F[(A{t)A)]. 

Where 

(i(t)i) = (V^|e+'^*/'"Me'*^*/''i|7/;„) (7) 



J{x{t),x{0)) ^e'^jy""'^' (9) 

is the Jacobian for the transformation of a volume 
element dx{0) to dx{t) along path x{t), dx{t) = 
J{x{Q),x{t))dx{Q). Its time evolution follows directly 
from the continuity equation with initial condition 
J(x(0),x(0)) — 1. Likewise, S is the action integral inte- 
grated over the quantum trajectory connecting a;(0) and 
x{t). Using this, we can rewrite Eq.[7|as an integral over 
starting points for Bohmian trajectories: 



(i(i)i) = / dx{0) / dy(0)A(x(i),2/(t))jy2jj;/2e-*^("(°)'"(*»/''V'*(a;(0))e*^-^(^(*')'^(*»/'^i^A(2;(0)) 



(10) 
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where = A\'4') a-nd we assume that A is a scalar 

quantum operator with matrix elements (a;|^|y) = 
A{x^y). In this expression, the forward-going quantum 
path along y{Q) to y{t) is subject to a quantum potential 
derived from the evolution of '0a (i) (with action inte- 
gral SA{y{Q),y{t)) and Jacobian Jy), where as the re- 
verse path along a;(0) to x{t) has a quantum potential 
derived from the evolution of i^{t) (with action integral 
S[x{Q),x{t)) and Jacobian J^.) Thus, even if x(0) = y{Q) 
the forward path {y{s)) will be very different than the re- 
verse path. We can interpret this in the following way. 
At time t — Q, we act on the initial wave function with 
operator A. If -0 is not an eigenstate of A, then we have 
a new wave function A\il)) — IV'a) which evolves out to 
time t where we act again with A and y{t) is instantly 
scattered to x{t) where it evolves back in time under the 
reversed evolution of V'(^)- The IVR correlation function 
then measures the probability amplitude that the final 
forward-backward state at x(0) was a member of the ini- 
tial state at y(0). This is illustrated schematically in 
Fig.d 

So far, no approximations have been made. The sim- 



plest correlation function is that the autocorrelation 
function of the wavefunction which monitors the overlap 
between the initial wavefunction ipo and the time-evolved 
wavefunction. 



c{t) = {mm)) 

= / dx{Q)ro{x{t)Mm)j^'^e''"'' 
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Here, ij^oix) is the initial wave function. This is directly 
analogous to the semi-classical IVR develo ped exten- 
sively over the past number of years. [H HI |33 
In the SC version, one estimates the quantum propagator 
via the well known Van Vleck expression ^25] . 



cit) « j dx{Q)ro{x(t))Ksc{x{t),x{Q)-t)i;o{x{Q)m) 
Here, Ksc{x{t),x{())\t) is the semi-classical propagator 



Ksc= 



cl. paths 
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A(x(t),y(t)) 




FIG. 1: Schematic representation of how correlation functions 
are constructed within the Bohm initial value representation 
of Ea llOl We start off with an initial wavefunction, ■(/ja = AiIj, 
sample an ensemble of starting points {y(0)} and evolve quan- 
tum trajectories out to time t. Simultaneously, one evolves 
a series of quantum trajectories backwards in time starting 
with ■(/)* as the initial wave function. The two trajectories are 
connected at time t via the matrix elements of A at the end 
points, A{x{t),y{t)). 



where the summation is over all classical paths connect- 
ing a:(0) to x(t) with classical action Sd, 



C is the matrix of the negative second variations of the 
classical action with the end points, 



C 



dx{0)dx{t) 



(15) 



and K gives the Morse index obtained by counting the 
number of conjugate points along each trajectory. What 
is appeaUng in the Bohm-IVR expression above is that 
one does not need to compute C or k along each tra- 
jectory. This is not surprising since we are dealing with 
a purely quantum mechanical object and C and k are 
semi-classical quantities. 

A closer connection between Eg. 1121 and Eq.^Jcan be 
drawn by writing C as 



C = - 



dpjt) 
dx{0) 



= ~J 



dx{tf 



(16) 
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(14) Thus, the SC-IVR expression can be written as 



Csc{t) 



dx{0) 



E 

cl. paths 



jl/2 



dx^t) 



1/2 



(17) 



Comparing Eq.^lto Eq.^]we first note the appearance 
of J^/^ in both expressions. Here, as in the case above, 
J is the Jacobian for the transformation of a volume ele- 
ment dx{0) to dx{t) along the classical path Xci(t). What 
is most striking, however, is that in the exact Bohm IVR 
expression one has a single trajectory connecting each ini- 
tial point to each final point where as in the semi-classical 
case, we have to consider every possible trajectory that 
starts at x{0). In the semi-classical case, this can take the 
form of a phase-space integration, as in the wi dely used 
Herman-Kluk expressionll IM EHH El Hfl |3l| . 

As mentioned above, one of the primary difficulties en- 
countered in developing numerical methods based upon 
the Bohm trajectory approach is in computing the quan- 
tum potential, Q, and its derivatives on a generally 
unstructured mesh of points evolving according to the 
Bohmian equations of motion. Moreover, the trajecto- 
ries themselves can be highly erratic and kinky requiring 
that one use very accurate numerical integration methods 
to achieve an accurate solution of the equations of mo- 
tion. |l7j | Consequently, obtaining an accurate estimate of 



the time-evolving wavefunction (via the synthesis equa- 
tion above) is quite difficult. On the other hand, one 
may be able to get by with a less robust and less ex- 
act calculation of the quantum trajectories if one is only 
interested in computing time-correlation functions. We 
next turn our attention towards a hierarchy approach for 
generating approximate quantum trajectories. 



III. DERIVATIVE PROPAGATION 

Recently, Trahan, Hughes, and Wyatt '2l'| introduced 
an interesting approach for computing quantum trajec- 
tories by simultaneously evolving the spatial derivatives 
of both C and S along a given trajectory. This notion 
is similar in spirit to the hydrodynamic moment expan- 
sions of the density matrixfj H 0, 0, 0| and follows 
directly by taking derivatives of dC / dt and dS/ dt with 
respect to x and transforming to the moving Lagrangian 
representation. Following this notion, the time-evolution 
of the n-th spatial derivatives of C and S at Xk{t) are 
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given by 




Our notation is that fl" = V"/ evaluated at Xk{t) and 
(") is the binomial coefficient. These equations, along 

with dxj/dt = S^^"^ /m are an infinite hierarchy of equa- 
tions for both C and S and their derivatives evaluated 
along the path Xj (t) . In this hierarchy, the equation for 
the n*^ order derivative of C involves up to the n + 2 and 
order derivative of S and up to the n + 1 order deriva- 
tive of C. Likewise, the equation of motion for the n*'* 
order derivative of S involves up to 2 higher order terms 
of C and one higher order term of S. The C^""* and S'["'' 
terms are in essence n*''-order cumulants of the wave- 
function evaluated at Xk{t). Note also that the hierarchy 
equations are local in that they do not couple different 
trajectories. Each trajectory can be computed as an in- 
dependent quantity and the task of computing multiple 
quantum trajectories can be easily distributed over mul- 
tiple machines. 

Imbedded in these equations of motion is the mech- 
anism by which complex structure emerges within the 
wave function since even simple wave functions with few 
cumulants will rapidly evolve into wave functions with 
high order cumulants. Thus, in order to close the equa- 
tions of motion, we need to introduce an artificial bound- 
ary condition which limits the extent of the hierarchy. 
The most straight-forward scheme is to simply require 
any term, cj"' or S'j"'' with n > riorder to be set to 
zero in the initial equations of motion. One can also 
introduce "absorbing boundary conditions" by damping 

higher-order terms, e.g. dSj"''' /dt = —XSj"\ 

To compute correlation functions, we first start off by 
sampling ipo on a uniform grid of DVR points such that 
■0071 = tpoixn{0))wl/'^ where wl/"^ is the gaussian quadra- 
ture weight associated with the n-th grid point. Thus, 
the correlation function can be evaluated by gaussian 
quadrature over a grid defined by the initial positions 
of the Bohm trajectories. For example, for the auto- 
correlation function of the wavefunction: 

Cit) = J2 WnroM))M^nmJ^/^e'^"/^ (20) 
n 



where J„ and Sn are the Jacobian and quantum action 
along Xn{s) which starts at Xn{0) and ends at Xn{t). 

Written in terms of the derivative terms above, the 
wave function correlation becomes simply 

N 

c{t) = 5]«;,e-<'(%<'W/''-^:(x,(O)0o(x,(O)) (21) 
i=i 

where cj"^ and 5'j"^are computed using the above hi- 
erarchy equations along a finite set of trajectories. As 
each trajectory is independent, the sum can evaluated by 
distributing the trajectory calculations amongst several 
nodes, each running a single quantum trajectory. This 
last expression along with its implementation comprise 
the central results of this paper. In the following section 
we examine a series of one dimensional test cases. 



IV. SAMPLE CALCULATIONS 

As discussed above, quantum interference in the wave 
function leads to very complex structure in the quan- 
tum potential. Even in the simple double-slit experiment 
example, the quantum potential due to the interference 
between the two gaussian components of the wave func- 
tion varies rapidly as the underlying quantum trajecto- 
ries avoid one another. 201] Consequently, in numerical 
applications of the quantum trajectory approach, nodes 
and interferences in the wave function have present con- 
siderable hurdles in the development of generalized multi- 
dimensional applications of this method. 

Here, we take three related test cases with increasing 
anharmonicity and test the ability of the Bohm-IVR ap- 
proach with derivative propagation to obtain accurate 
correlation functions. In each case we consider an elec- 
tron (m = 1) in one-dimensional potentials, chosen such 
that the curvature at the bottom is identical in each case. 

Vho = x'/2-l (22) 
Vqo = a;V2 + O.Ola;'' - 1 (23) 
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FIG. 2: Potential functions used in comparative examples. 
(Key: — — Harmonic (Vho); ■ ■ ■ — Quartic, (Vqo); — ■ 
— =Gaussian well (Vg).) 

Vg = -e-"'/2 (24) 

Vho is the harmonic well, Vqo is harmonic well per- 
turbed by a small quartic term, and Vg is an inverted 
gaussian well. In each case, the initial state '4^{x) — 
exp(c + is/h) is gaussian with c and s given by 

c{x) = i log(/3/7r) - (3(x - xof (25) 
s{x) ^ (26) 

with Xo = 1 and j3 = 1. (Note, atomic units used 
throughout.) A plot of each potential along with the 
location of the lowest energy bound states is show in 
Fig.H 



A. Harmonic Oscillator 

Harmonic systems are particularly useful as test cases 
since the an initial gaussian wave function retains its 
gaussian form for all time. Furthermore, any quantity 
of interest can be computed analytically thereby provid- 
ing a useful benchmark for the method. In Fig. |3 we 
show the quantum trajectories computed via the deriva- 
tive propagation scheme outlined above. The trajectories 
display the coherent oscillation one expects for Bohm tra- 
jectories in this system. This is not surprising since a low 
order truncation of the derivative propagation equations 
is perfectly valid for this case. These agree perfectly with 
the analytical results for this system and we can inte- 
gration the derivative propagation equations out to very 
long times without significant loss of accuracy. The initial 



x(t) 




FIG. 3: Quantum trajectories for harmonic oscillator com- 
puted using derivative propagation approach. 

points were chosen as evenly spaced gaussian quadrature 
points. 

The solid curve in Fig. 0] shows the time correlation 
function of the wavefunction and its Fourier transform. 
Here, we performed the calculation out io t = 20n 
corresponding to 10 complete oscillation cycles. Both 
C{t) and its spectral representation C(u;) agree per- 
fectly with analytical results with C{u!) accurately map- 
ping out the harmonic progression of the energy lev- 
els in this well ojn — w(n + 1/2) as well as the ex- 
pansion coefficients for the projection of the initial 
wavefunction on to the eigenstates of the well, c„ = 
KnlV^o)!^- For the case at hand, the first 4 coefficients 
are c = {0.606531,0.303265,0.0758163,0.0126361} and 
agree nicely with the roughly 8:4:1 ratio of the first three 
peaks in C(w). 



B. Quartic Perturbed Harmonic 

In this case, we apply a slight anharmonic squeeze to 
the harmonic well. In an anharmonic system, an ini- 
tial gaussian wave function retains its gaussian form for 
short period of time before bifurcating into a multitude of 
components due to interferences. Such features are most 
likely to occur where the underlying classical trajectory 
crosses a conjugate or focal point. Such crossings are 
strictly forbidden for quantum Bohmian trajectories and 
lead to rapidly varying quantum forces over a short spa- 
tial region. Consequently, interference effects in anhar- 
monic systems has proven to be somewhat of a bugbear 
for implementing quantum trajectory approaches based 
upon fitting the quantum potential (and its derivatives) 
to a simple polynomial over a set of quantum trajectories. 

The derivative propagation trajectories also retain 
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FIG. 4: Correlation function and its spectral representation 
for a gaussian wavepacket in various potentials computed us- 
ing Bohm-IVR/derivative propagation approach. (Key; — = 
Harmonic (Vho)', ■ ■ ■ — Quartic, (Vqo); — ■ — =Gaussian well 
(Vg).) 



their coherences for a short period of time. However, 
after only one or two oscillations, crossings become more 
frequent and the ensemble of trajectories degrades into 
an incoherent mess. Had we performed this with the 
moving-weighted least squares (MWLS) or other meth- 
ods which extracts the quantum potential from an en- 
semble of trajectories, the results would have been cat- 
aclysmic. Our experience has indicated that even small 
errors in constructing the quantum force from the en- 
semble can lead to a rapid amplification of error in the 
propagation. Remarkably, even though the ensemble has 
apparently lost coherence amongst its members, the indi- 
vidual trajectories exhibit non-classical deflections, and 
quasi-periodic behavior one expects from truly quantum 
trajectories. 




FIG. 5: Quantum trajectories for inverted gaussian well ex- 
ample. 

Turning towards the autocorrelation and its Fourier 
transform. C{t) shows some recurrences at i « 27rn as 
in the harmonic case, but these rapidly die out after 5-6 
beats. This limits our ability to resolve the eigenenergies 
in the Fourier transform. Nonetheless, at least 4 peaks 
can be distinguished. 



C. Inverted Gaussian Potential 

In this final example, Vq is an inverted gaussian and 
supports only two eigenstates (at E = —0.59386 and 
E = —0.0356576). In this example we used 30 trajecto- 
ries since the majority of the trajectories near the leading 
edges of the distribution leave the potential region al- 
most immediately. In Fig. [3 we show 10 trajectories that 
started in the center of the initial wavepacket. As time 
progresses, many of these eventually escape the potential 
well leaving only one of the original 30 trajectories in the 
well. 

The long-time accuracy of the ensemble is not so im- 
portant since we are ultimately interested in computing 
C{t). This we show in Fig. 0] along with its Fourier trans- 
form. One can clearly identify 6 recursion peaks spaced 
roughly every 0.35 ps. These correspond to the nearly 
harmonic motion of the trajectories closest to the center 
of the ensemble. 



V. DISCUSSION 

For low dimensional systems, there are arsenals of 
robust and computationally efficient methodologies for 
computing exact quantum dynamics for both time de- 
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pendent and time-independent systems. However, these 
methods are severely hmited by dimensionaUty since the 
computational effort and storage required to perform 
an exact quantum calculation increases almost exponen- 
tially with dimensionality. Thus, approximate methods, 
such as those introduced here, provide a viable path way 
towards higher dimensions with increasing complexity. 
While we have focused entirely upon one-dimensional ex- 
amples, the Bohm-IVR and derivative propagation equa- 
tions are exact and can be easily implemented in a multi- 
dimensional setting. 

What do the individual trajectories mean? Are 
the trajectories computed via the derivative propagation 
methods really Bohm trajectories? The answer, I believe 
is a definite yes and no. Yes, they are Bohm trajectories, 
but not for a initial single gaussian wavepacket evolv- 
ing in a given potential. They correspond to trajectories 
arising from a diverging set of wave functions where each 
wave function component is constrained to have a par- 
ticular functional form dictated by truncating the hier- 
archy. In essence, by truncating the hierarchy, one intro- 
duces dephasing and decoherence into the dynamics of a 
quantum system. This is evident in the time-correlation 
function for the particle in the harmonic -I- quartic po- 
tential. The time correlation function for a fully coherent 
wavepacket should show stronger recurrences and even- 
tually be able to resolve the eigenenergies of the system 



upon Fourier transform. One evidence for this is that for 
bound systems other than the harmonic oscillator, the 
autocorrelations we compute using the derivative propa- 
gation decay to zero after a few recurrences, and never 
recover coherence. 

For exact quantum mechanics, this implicit loss of co- 
herence at long-times clearly undesirable. However, for 
condensed phase systems, especially in cases where we 
have a low dimensional quantum subsystem in contact 
with a high-dimensional "bath" of oscillators, reduction 
of the density matrix of the subsystem to a gaussian form 
in the coherences occurs spontaneously due to decoher- 
ence. dill El. Consequently, so long as the "physical" 
decoherence time is smaller than the "truncation induced 
" decoherence time, time-correlation functions computed 
using the Bohm-IVR approach we describe herein arc ex- 
pected to quite accurate. 
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